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density functional theory in real time: application to polycyclic aromatic 

hydrocarbons 

Miguel A. L. Marques, 1 ' 2 'B Alberto Castro, 3 ' 2 Giuliano Malloci, 4 Giacomo Mulas, 4 and Silvana Botti 5 - 1 ' 2 

1 Centro de Fisica Computational, Departamento de Fisica, Universidade de Coimbra, Coimbra, Portugal 

2 European Theoretical Spectroscopy Facility 
3 Institut fur Theoretische Physik, Fachbereich Physik der Freie Universitat Berlin, Arnimallee 14, D-1J.195 Berlin, Germany 

^INAF Osservatorio Astronomico di Cagliari, Astrochemistry Group, 
Strada n.54, Loc. Poggio dei Pini, 109012 Capoterra (CA), Italy 
5 Laboratoire des Solides Irradies, CNRS-CEA-Ecole Poly technique, Palaiseau, France 

(Dated: February 6, 2008) 

The van der Waals dispersion coefficients of a set of polycyclic aromatic hydrocarbons, ranging 
in size from the single-cycle benzene to circumovalene (C66H20), are calculated with a real-time 
propagation approach to time-dependent density functional theory (TDDFT). In the non-retarded 
regime, the Casimir-Polder integral is employed to obtain C§, once the dynamic polarizabilities have 
been computed at imaginary frequencies with TDDFT. On the other hand, the numerical coefficient 
that characterizes the fully retarded regime is obtained from the static polarizabilities. This ab 
initio strategy has favorable scaling with the size of the system - as demonstrated by the size of the 
reported molecules - and can be easily extended to obtain higher order van der Waals coefficients. 
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I. INTRODUCTION 

When two molecules are far apart, they still inter- 
act through long range electromagnetic forces, named 
after J. D. van der Waals X These are crucial to under- 
stand a wide range of phenomena: the chemistry of rare 
gases. 2 protein folding and dynamics, 3 the machinery of 
neurotransmitters^ and the molecular chemistry in the 
interstellar medium- are just some examples. In partic- 
ular, in astrochemistry van der Waals forces are essen- 
tial to describe neutral- neutral reactions, which are not 
yet well accounted for in rate-reaction models. These 
intermolecular forces are mainly originated by electric 
multipole-multipole interactions. Even if the static mul- 
tipoles of the molecules are null (as the dipole of a single 
spherical atom in the ground state), they still contribute 
to the interaction due to their spontaneous oscillations. 
These contributions of dynamical origin, that cannot 
be explained by Electrostatics, are the London disper- 
sion forces £ These are attractive interactions, and arise 
from the mutual polarization of the two electronic clouds. 
They dominate the long range dynamics of the molecules, 
and are also present at short distances - although masked 
by other stronger (ionic or covalent) forces. 

The calculation of these dispersion forces can be a 
challenging problem^ Three regimes have to be distin- 
guished, according to the distance that separates the in- 
teracting molecules: 

(i) short distances, such that there is non-negligible 
overlap of the electronic clouds of the two molecules. 
This is the most difficult situation, since it requires, in 
principle, a supermolecule calculation, i.e. the treat- 
ment of the two molecules combined together as a sin- 
gle entity. The valid approaches in this overlapping 
regime (and therefore also valid in the non-overlapping 
regime) are sometimes called seamless van der Waals 



techniques. Numerous possibilities exist, although none 
of them entirely satisfactory: full configuration interac- 
tion for very small molecules,— M0ller-Plesset perturba- 
tion theory^ or Monte Carlo methods* 1 ^ Attempts to 
use ground state density functional theory (DFT) are 
specially challenging^ In the realm of time-dependent 
density- functional theory (TDDFT), the recent work of 
Dobson 11 summarizes the current research, largely based 
on the adiabatic connection / fluctuation dissipation the- 
orem. Dispersive forces can also be added to standard 
DFT through empirical correction terms. These, how- 
ever, require the previous knowledge of the Cg van der 
Waals coefficients (see below), often roughly estimated 
from the atomic Cq's .— In the present work we will not 
consider this regime. 

(ii) long distances, such that we can neglect the over- 
lap. In this case, the electrons belonging to different 
molecules are distinguishable, and one can isolate the 
Coulomb operator that corresponds to interactions be- 
tween electrons of different molecules, and apply (second 
order) perturbation theory for this operator J£ The first 
term in the perturbative expansion of the interaction en- 
ergy decays as — Cg/i? 6 , where R is the intermolecular 
distance. 

(iii) very long distances, such that retardation effects 
become important ^ This means that the time it takes 
for the photons that mediate the electromagnetic inter- 
action to travel between the two molecules is not negligi- 
ble. Retardation is described by a correction factor that 
is equal to unit for small distances and proportional to 
1/R for large distances. 

The knowledge of the various static and dynamic mul- 
tipole polarizabilities at imaginary frequencies suffices to 
compute the van der Waals interaction in the regimes 
(ii) and (iii). These response functions can be calculated 
within a variety of quantum chemical methods, typically 



2 



through the evaluation of a wide range of excitation en- 
ergies and associated transition matrix elements. Alter- 
natively, at a lower computational cost, we can use time- 
dependent density- functional theory (TDDFT)^^ This 
theory is very successful in predicting optical spectra (i.e., 
essentially the imaginary part of the dynamic polarizabil- 
ity at real frequencies) for finite systems, and has been 
used to study a wealth of molecules and clusters M> The 
first calculation of Cq coefficients using TDDFT was per- 
formed in 1995 by van Gisbergen et al for a variety of 
small molecules »i£ 

TDDFT in the perturbative regime can be formulated 
in a variety of flavors*^ With the purpose of calculat- 
ing Cq coefficients, one can find several approaches. Van 
Gisbergen et al» opted for a self-consistent calculation 
of the response function. The matrix formulation of lin- 
ear response theory 2 ^ can also be employed, yielding os- 
cillator strengths and excitation energies, sufficient in- 
gredients for the computation of the dynamic polariz- 
ability. A slightly different approach is the polariza- 
tion propagator technique^ which can be constructed 
on top of TDDFT (or on top of other electronic struc- 
ture theories, e.g. time-dependent Hartree-Fock, TDHF). 
The linear polarization propagator was demonstrated to 
work for complex frequencies^ opening the possibility of 
computing Cq coefficients^ With this technique, Jiem- 
chooroj and collaborators computed Cq coefficients of no- 
ble gases atoms j^ 3 - n-alkanes^ 3 - (n < 7), polyacenes^ the 
Ceo molecule^ and sodium clusters up to 20 atoms^ 
For large molecules, Banerjec and Harbola have proposed 
the use of orbital free TDDFT, 26 providing satisfactory 
results for large sodium clusters. 

In this work, we propose an alternative scheme, based 
on the explicit propagation of the time-dependent Kohn- 
Sham equations ! 27 ' 28 This approach has already proved 
its usefulness to calculate the dynamic polarizability of 
large systems (see, e.g., Reff29l). The explicit time- 
dependent picture, where one is confronted with an ini- 
tial value problem, can be preferable over the the time- 
independent picture - where the mathematical problem 
is diagonalization -, especially regarding scalability with 
the size of the s yste m. For a discussion in the field of 
TDDFT, see Ref[jJ. 

We have chosen an extensive set of polycyclic aromatic 
hydrocarbons (PAHs) in order to validate the methodol- 
ogy. PAHs, a large class of conjugated 7r-electron sys- 
tems, are of great importance in many areas, among 
which combustion and environmental chemistry, mate- 
rials science, and astrochemistry. PAHs are found in 
carbonaceous meteorites, in interplanetary dust parti- 
cles, and are thought to be the most abundant molecular 
species in space after H2 and CO, playing a crucial role in 
the energy and ionization balance of interstellar matter 
in galaxies. Motivated by this astrochemical relevance, 
some of us have recently performed an extensive study of 
these systems! 30 ' 31 This research is collected in a thor- 
ough compendium of molecular properties of PAHs. 31 
The computation of van der Waals constants for pairs 



of PAHs, as a first approach to the analysis of their in- 
tcrmolecular properties, is therefore a natural extension 
of this work. In fact, van der Waals parameters are a 
crucial ingredient to model condensation and evapora- 
tion of PAH clusters, which are another important player 
in the physics/chemistry of the interstellar medium. In 
current works, people usually employ empirical van der 
Waals parameters for the relatively long-range part of 
the interaction in molecular dynamics simulations, to- 
gether with some tight-binding approximation for the 
short-range parti 3 ^ 

One previous calculation of dispersion coefficients of 
PAHs, based on the complex linear polarization propa- 
gator method, was reported in Reff24l. although limited 
to benzene, naphthalene, anthracene, and naphtacene. 
This fact adds a further motivation to our study, since it 
permits to compare the results of the two schemes. 

II. METHODOLOGY 

When both the wavefunction overlap and the retarda- 
tion effects are almost null [situation (ii)] , the application 
of second order perturbation theory leads to an expan- 
sion of the interaction energy with respect to the inverse 
of the intermolecular distance (l/R): 

n— 6 

The coefficients C„ are usually called Hamaker 
constants^ 3 - The leading non-null term Cq is due to the 
dynamic dipole-dipole polarizability. The odd terms un- 
til n — 11 are also null for spherical molecules (or if an 
average over relative orientations is taken) , if we neglect 
retardation effects. In principle, the coefficients also de- 
pend on the relative orientation of the molecules. 

The Cq B dispersion Hamaker constant for a pair of 
molecules A and B, averaged over all possible relative 
orientations, is related to the dipolc molecular polariz- 
ability through the Casimir-Polder relation (atomic units 
will be used hereafter): 

o /-oo 

C£ B = - / du a< A >(iu) a^iiu) , (2) 

T JQ 

where a^ x \\u) is the average of the dipole polarizability 
tensor of molecule X, o^ x \ evaluated at the complex 
frequency iu: 

a W(iu) = -Ti[a (x \m)] . (3) 
3 

It should be noted that: (i) If we fix the relative orienta- 
tion of the molecules, the orientation dependent Hamaker 
constant can also be calculated from the a tensors, con- 
sidering the appropriate linear combination of their com- 
ponents; (ii) higher order Hamaker constants, useful for 
shorter distances, can be obtained through the use of 
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analogous formulae involving higher order multipole po- 
larizability tensors. The expressions accounting for these 
two generalizations can be found, for example, in RefflH. 
The calculations presented below are solely concerned 
with Eq. ([2]) ; however we wish to stress that the method- 
ology trivially yields full polarizability tensors of arbi- 
trary order, hence providing the possibility to tackle 
those general cases, with almost no extra computational 
cost. 

On the other hand, when the distance increases and re- 
tardation effects become important, the van der Waals in- 
teraction depends solely on the static polarizability ; 14 i 34 
decaying as AE(R) = —K AB /R 7 , where the constant is 



rAB = 23c (a )(0) q(b) 
8ir z 



(4) 



where c is the velocity of light in vacuum. 

The TDDFT time propagation approach to obtain 
the dynamic polarizability components can be summa- 
rized as follows: let us assume a weak electrical (spin- 
independent) dipole perturbation 
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r (r,co) = -XjK(tL>) . 



(5) 



This defines an electrical perturbation polarized in the 
direction j: 5E(a>) = K(ui)e.j. The response of the system 
dipole moment in the i direction 

8{Xi)(u) = J2[ d3rx * Sn -( r ^) (6) 
is then given by: 
8(£i)(u) = ~k{lo)J2 Jd 3 r fd 3 r'xi w(r,r' » x) , 

GO 1 

(7) 

where x<j<j' is the response function of the system: 

Sn a {r,w) = J d 3 rxaa'(r,r' ,uj)Sv cxt ,a'(r' ,uj) ■ (8) 

The dynamic dipole polarizability a,^{tJ) is the quotient 
of the induced dipole moment in the direction i with the 
applied external electrical field in the direction j, which 
yields: 

Oijiu) = -^2 Jd?rJ dV x { x*a> (r, r', w) x'j . (9) 

The dynamic polarizability elements may then be ar- 
ranged to form a second-rank symmetric tensor, a(u>). 
We consider now a sudden external perturbation at t = 
(delta function in time, which means k(u>) — ft, equal for 
all frequencies), applied along a given polarization direc- 
tion, say ij. By propagating the time-dependent Kohn- 
Sham equations,— we obtain 5n(r,t). The polarizability 
element atij(ui) may then be calculated easily via: 
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d r Xi 5n(r, oj) 



In order to obtain values of a at imaginary frequen- 
cies one only has to substitute to by iu. This computa- 
tional framework has been implemented in the octopus 
code, described in 

RefH| 

to which we refer the reader 
for technical details. The octopus code performs the 
calculation on a real space mesh, which reduces the con- 
vergence problem to two parameters: the grid spacing 
(we used 0.3 A in this case) and the size of the simula- 
tion box. The ion-electron interaction was modelled with 
norm-conserving Troullier-Martins pseudopotentials^i 

The geometries of the PAHs were obtained by means 
of DFT as described in Refl3ll For the optimiza- 
tion of the geometries, we chose the hybrid B3LYP 38 
as the approximation to the exchange and correlation 
functional. For the TDDFT propagations, however, we 
used the adiabatic local density approximation (LDA), 
which has proved to yield reliable results for conjugated 
molecules^ the use of more sophisticated functionals is 
possible, but does not change the results for the dynamic 
polarizability significantly^ 



III. RESULTS 

We have computed the Cq B Hamaker constants for 
all possible pairs {A, B} of a set of 41 PAHs. The homo- 
molecular Hamaker constants C§ A for all the PAHs stud- 
ied are reported in Tables Q] and [TTJ the full set can be 
consulted in the database described in Refl3~0. The ta- 
bles also display the average static dipole polarizability 
ct(0), the effective London frequency uj\ (see below) and 
the retarded van der Walls coefficient K . 

It is very difficult to extract Hamaker constants exper- 
imentally; to our knowledge, there are no experimental 
results reported for any PAH to compare to. For ben- 
zene, however, Kumar and Meith 4 ^ reported a value of 
1723 a. u., by making use of the dipole oscillator strength 
distributions (DOSDs), which are constructed from ex- 
perimental dipole oscillator strengths and molar refrac- 
tivity data. This is in good agreement with our computed 
value of 1763 a. u. 

We have also displayed for comparison the numbers 
reported in ReflM obtained by means of the complex 
linear polarization propagator scheme. This scheme was 
constructed on top of TDDFT, although making use of 
the B3LYP functional. These methodological differences, 
and further numerical details can explain the very small 
differences, in all cases below 2% - see the results in 
the Table |T] for benzene, naphthalene, anthracene and 
tetracene (also known as naphthacene) . 

In the so-called London approximation, the polariz- 
ability at imaginary frequencies is modelled with the help 
of only two parameters: the static polarizability, and one 
effective frequency lo\. 



a(iu) 



a(0) 



(11) 



-- jdt j d i rx i 5n(r,t)e 



(10) 



i + W^i) 2 ■ 

Upon substitution on the Casimir-Polder integral, this 
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yields for a homo- molecular Hamaker constant: 

a = ^V(o). 



(12) 



With the knowledge of Cg and a(0), one can obtain u>\ 
from Eq. [12] These effective frequencies are also re- 
ported in Tables U and HH They are roughly decreas- 
ing with the size of the PAH, from 0.482 Ha for ben- 
zene to the 0.156 Ha of pentarylene. The decrease is, 
however, strongly irregular. As it has been pointed out 
before,— it can be related to the ionization potential of 
the molecules; this is demonstrated in Fig. [1] where we 
have plotted the ionization potentials of the PAHs (cal- 
culated at the DFT/B3LYP level of theory), versus the 
effective frequency u>%. The data points approximately 
accumulate around a straight line, proving the correla- 
tion. 

Equation 1121 also gives us a hint on the dependency 
of Ce with the size of the molecule: it is proportional 
to the square of the polarizability (the product of polar- 
izabilities, if the molecules are different), which in turn 
typically grows with the volume, and therefore, with the 
number of atoms. Consequently, one should expect a lin- 
ear dependency of C AB with respect to the product of the 
number of atoms, Na x Nb- This is indeed confirmed in 
Fig.O Some cases, however, deviate from a straight line. 
These cases correspond to strongly anisotropic PAHs 
(with three very different axes). This is captured by the 
dipolc anisotropy: 



Aa 2 
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(13) 

Figure [3] shows how the PAHs whose Hamaker constant 
deviates strongly from the general trend are those whose 
polarizability anisotropy is also stronger: we have over- 
laid the values of Cq for all PAHs (divided by the square 
of the number of Carbon atoms N 2 ), with the values Aa 2 
(also divided by N 2 ). One can see how the two datasets 
are correlated - specially in the right side of the graph, 
which corresponds to the larger PAHs. 

The crossover between the non-retarded and retarded 
regimes is given by the length scale A = 2ttc/0j, where 
u> is a characteristic frequency of the electronic spec- 
trum of the molecule. For R/X > 10 we enter the fully 
retarded regime, while for R/X < 0.1 we can still use 
Eq. ([T]). Using for Q the values of uj± obtained through 
the London approximation, we reach values in the range 
of A ~ 0.1 /im (for benzene) to A ~ 0.3 /im (for pentary- 
lene). 



It is interesting to notice that in the fully retarded 
regime we can write K AA as a function solely of the 
Hamaker C AA coefficient and the effective London fre- 
quency lo\. Combining equations ^ and (|12p we arrive 
at 



k aa _ 23c C 6 



6tt 2 



(14) 



The values of homomolecular coefficients K AA for the 
PAHs studied in this work can be found in Tables HI and [Hi 
IV. CONCLUSIONS 

In this Article we presented a method to calculate the 
van der Waals coefficients of molecular systems using a 
time-propagation scheme within TDDFT. As an exam- 
ple, we have applied this method to the family of poly- 
cyclic aromatic hydrocarbons. Our results are in excel- 
lent agreement with available theoretical and experimen- 
tal data, and fully validate our approach. Values of Cq 
scale approximately with Na X Nb , where Nx is the num- 
ber of atoms in the molecule X. The strongest deviations 
from this law are for highly anisotropic structures. 

This scheme has several non-trivial advantages: (i) It 
scales with N 2 , where _/V a is the number of atoms in the 
molecule, (ii) The time-propagation yields the polariz- 
ability in real time. From this quantity it is then im- 
mediate to obtain the real and imaginary parts of a at 
real and imaginary frequency. Therefore, the optical ab- 
sorption spectrum, static polarizabilities, Cg coefficients, 
etc. are calculated in one shot, (iii) This scheme is easily 
generalized to higher order coefficients and arbitrary ge- 
ometries. We therefore expect it to be useful in the study 
of large systems, like clusters or molecules with biological 
interest. 
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TABLE I: Average static polarizability a(0), Ce Hamaker 
constant, effective frequency uj\ — (4/3)Ce/a 2 (0), and the 
coefficient K of the retarded van der Waals interaction for 
the selected PAHs. All quantities are given in atomic units. 



Molecule 


a(0) 


C 6 


Wl 


K 






xlO -3 




x 10 -6 


benzene (CeHe) 


70.5 


1.76 


0.473 


0.198 


(RefUl TDDFT/B3LYP) 


70.0 


1.77 


0.482 


0.195 


(ReflU DOSD) 


67.8 


1.72 


- 


- 


azulene (CioHs) 


133 


5.15 


0.390 


0.703 


naphthalene (CioHs) 


123 


4.79 


0.422 


0.605 


(Ref[2l; TDDFT/B3LYP) 


122 


4.87 


0.439 


0.594 


acenaphthene (C12H10) 


143 


6.76 


0.439 


0.820 


biphenylene (C12H8) 


152 


6.91 


0.400 


0.920 


acenaphthylene (C^Hg) 


145 


6.57 


0.417 


0.838 


fluorene margin (C13H10) 


159 


7.97 


0.422 


1.00 


phenanthrene (C14H10) 


182 


9.36 


0.379 


1.32 


anthracene (C14H10) 


189 


9.92 


0.372 


1.42 


(ReflU; TDDFT/B3LYP) 


185 


10.0 


0.392 


1.37 


pyrene (Ci 6 Hi ) 


205 


12.0 


0.380 


1.68 


tetracene (C18H12) 


264 


17.5 


0.336 


2.78 


(ReflU; TDDFT/B3LYP) 


259 


17.5 


0.349 


2.68 


triphenylene (C18H12) 


231 


15.7 


0.390 


2.14 


benzo[a]anthracene (C18H12) 


246 


16.6 


0.364 


2.42 


chrysene (Ci 8 Hi 2 ) 


239 


15.9 


0.373 


2.27 


benzo[e]pyrene (C20H12) 


260 


18.9 


0.375 


2.69 


perylene (C20H12) 


262 


18.6 


0.361 


2.74 


benzo[a]pyrene (C20H12) 


277 


19.8 


0.345 


3.06 


corannulene (C20H12) 


244 


17.5 


0.390 


2.38 
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TABLE II: Continuation of Table [fl 



Molecule 


a(0) 


c 6 




K 






xlO" 3 




xlO" 6 


anthanthrene (C22H12) 


304 


23.5 


0.340 


3.69 


benzo[g,h,i]perylene (C22H12) 


282 


22.2 


0.373 


3.17 


pentacene (C22H14) 


353 


28.1 


0.301 


4.98 


dibenzo[b,def]chrysene (C24H14) 


355 


30.4 


0.321 


5.04 


coronene (C24H12) 


318 


27.0 


0.357 


4.03 


hexacene (C26H16) 


454 


42.1 


0.272 


8.23 


dibenzo[cd,lm]perylene (C26H14) 


384 


34.8 


0.314 


5.89 


bisanthene (C28H14) 


402 


37.6 


0.310 


6.45 


benzo[a]coronene (C28H14) 


386 


37.9 


0.339 


5.96 


dibenzo[bc,ef]coronene (C30H14) 


431 


44.0 


0.316 


7.42 


dibenzo[bc,kl]coronene (C30H14) 


459 


46.3 


0.293 


8.40 


terrylene (C30H16) 


484 


47.8 


0.272 


9.36 


ovalene (C32H14) 


453 


49.6 


0.323 


8.18 


tetrabenzocoronene (C36H16) 


565 


65.4 


0.273 


12.8 


circumbiphenyl (C38H16) 


538 


69.8 


0.321 


11.6 


circumanthracene (C40H16) 


612 


82.5 


0.293 


15.0 


quaterrylene (C40H20) 


799 


97.0 


0.203 


25.5 


circumpyrene (C42H16) 


631 


89.0 


0.298 


15.9 


hexabenzocoronene (C42H18) 


590 


85.9 


0.330 


13.9 


dicoronylene (C48H20) 


770 


122 


0.274 


23.6 


pentarylene (C50H24) 


1196 


168 


0.156 


57.1 


circumcoronene (C54H18) 


840 


150 


0.284 


28.1 


circumovalene (C66H 2 o) 


1099 


236 


0.260 


48.2 
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